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I.  INTRODUCTION  / 

The  in  flight  characteristics  of  an  artillery  shell  are  of 
major  importance  to  the  shell  designer,  ballisticion  and  ultimately  the 
artillery  field  commander  whose  mission  is  deployment  of  timely  and 
accurate  fire  power.  The  aerodynamic  properties  of  artillery  shell, 
such  as  pitching  moment,  Magnus  moment  and  drag  are  critical  to  the 
stability  of  shell  which  in  turn  significantly  affects  accuracy  and 
time  of  flight.  The  capability  to  determine  the  aerodynamics  of  shell 
is  required  over  a  wide  range  of  flight  regimes  since,  depending  on 
initial  launch  velocities,  artillery  shell  are  subject  to  subsonic, 
transonic,  and  supersonic  flight.  Projectile  aerodynamics  over  these 
various  flight  regimes  have  been  found,  in  some  cases,  to  change  by  an 
order  of  magnitude.  The  solution  techniques  utilized  must  therefore 
be  capable  of  computing  these  changes. 

concentrated  theoretical  and  experimental  research  program 
has  been  ongoing  at  BRL  in  order  to  develop  the  predictive  capabilities 
required  for  determining  projectile  aerodynamics.  Supersonic 
computations  using  combined  inviscid  flow  field  and  boundary  layer 
techniques  have  been  developed*^  Sturek  -1  ^ ,  et  -al  ?,  j*for  cone- cylinder 
and  ogive-cylinder  configurations!  Recent  results  have  been  obtained 
in  supersonic  flow  over  a  typical  boattailed  projectnes-by— Sehi-ff-and^ 
Sfcurdc'^'tising  modern  computational  techniques  for  solving  the  thin- 
layer  Navier-Stokes  equations. 

V  * 


0  1980 


80 


NIETUBICZ  $  KAYSER 


.Inviscid  transonic  computational  results  have  been  obtained 
by  Reklis''3  ,  et  al.,  for  a  secant-ogive-cylinder-boattail  shape.  The 
inviscid  techniques  give  good  results  for  pitch  plane  aerodynamic 
coefficients  at  small  angle  of  attack.  However,  this  technique  lacks 
the  ability  to  include  viscous  effects  at  transonic  speeds  and  thus 
compute  the  Magnus  effect  on  projectiles.  Techniques  which  have  been 
applied  in  supersonic  flow  for  combining  inviscid  and  boundary  layer 
methods  have  not  been  fully  established  for  transonic  flow.  These 
methods,  which  have  shown  good  results  for  ogive-cylinders  in 
supersonic  flow  at  low  angle  of  attack,  are  not  accurate  in  modeling 
the  severe  flow  expansion  in  the  vicinity  of  surface  discontinuities 
such  as  those  that  occur  at  the  cylinder-boattail  junction. 

The  solution  of  the  thin- layer  Navier-Stokes  equations,  which 
allows  for  the  simultaneous  computation  of  the  inviscid  and  viscous 
regions,  eliminates  the  need  for  matching  two  different  solutions. 
Additionally,  since  all  three  momentum  equations  are  retained,  the 
ability  to  compute  in  regions  of  separated  flow  is  achieved.  This 
paper  describes  the  governing  three-dimensional  thin-layer  Navier-Stokes 
equations  used  for  computing  flow  over  projectile  shapes  at  angle  of 
attack.  Secondly,  the  generalized  axisymmetric  formulation  used  for 
computations  at  a  =  0°  will  be  described.  A  description  of  the 
numerical  algorithm  and  results  will  follow.  Experimental  and 
computational  results  will  be  presented  for  a  secant-ogive-cylinder- 
boattail  projectile  shape  at  a  =  0°  and  a  =  2°.  Computational  results 
will  also  be  shown  for  a  ring  airfoil  shape  thus  demonstrating  the 
general  geometry  capability  of  the  present  numerical  scheme. 

II.  GOVERNING  EQUATIONS 

The  general  three-dimensional  thin- layer  Navier-Stokes 
equations,  used  for  all  cases  where  a  =  0°,  are  described  in  Section  I  la. 
The  thin- layer  generalized  axisymmetric  equations,  which  are  a  special 
case  of  the  3-D  equations,  are  described  in  Section  lib. 

a.  Three-Dimensional  Equations 


The  transformed  three-dimensional  thin- layer  Navier-Stokes 

equations 
written  as 


in. non-dimensional  and  strong  conservation  law  form  arc 
v4  ) 


3  q  +  3_E  +  3  1;  +  3  G  =  Re'1 3  S 
C  n  C  C 


tP> 


where  general  coordinate  transformations 
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£  =  £(x,y,z,t)  -  longitudinal  coordinate 


n  =  n(x,y,z,t)  -  circumferential  coordinate 
?  =  C(x,y,z,t)  -  near  Normal  coordinate 
t  =  t  -  time 


are  used  and 
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The  velocities 


U  =  £.+£u  +  £v  +  £w 
t  x  y  z 

V  =  n+  nu  +  nv  +  nw 

t  x  y  z 

W  =  £  +  Cu+£v+£w 
t  x  y  z 

represent  the  contravarient  velocity  components. 


(2) 


The  Cartesian  velocity  components  (u,v,w)  are  retained  as  the 
dependent  variables  and  are  nondimensionalized  with  respect  to  a  (the 
free  stream  speed  of  sound).  The  local  pressure  is  determined  using 
the  relation 


P  =  (Y  -  1)  (e  -  .5p  (u2  +  v2  +  w2))  (3) 

where  y  is  the  ratio  of  specific  heats,  density  (p),  is  referenced  to 
p  and  total  energy  (e)  to  p  a  2.  The  additional  parameters  are  (k) 
tKe  coefficient  of  thermal  conductivity,  (y)  the  dynamic  viscosity, 

(Re)  the  Reynolds  number,  (Pr)  the  Prandtl  number,  and  (X)  which 
through  the  Stokes  hypothesis  is  (-2/3)y. 

The  metric  terms  £  ,  n  and  £  are  formed  from  the  derivatives 
x’  x  x 

x  ,  y  ,  z  ,  etc.,  and  together  with  the  Jacobian  of  the  transformation 

S  S> 

allow  for  computations  to  be  performed  for  variable  body  geometries. 

The  "thin- layer"  approximation used  here  requires  that 
all  body  surfaces  be  mapped  onto  £  =  constant  planes  and  that  Re  >>  1. 
Essentially,  all  the  viscous  terms  in  the  coordinate  directions  (here 
taken  as  £  and  n)  along  the  body  surface  are  neglected  while  terms  in 
the  £  or  the  near  normal  direction  to  the  body  are  retained.  This 
approximation  is  used  because,  due  to  computer  speed  and  storage 
limitations,  fine  grid  spacing  can  only  be  provided  in  one  coordinate 
direction  (usually  taken  as  the  near  normal  direction)  and  the  grid 
spacing  available  in  the  other  two  directions  is  usually  too  coarse  to 
resolve  the  viscous  terms.  For  the  type  of  problems  currently  under 
investigation,  i.e.,  projectiles  at  low  angles  of  attack,  with  no  strong 
cross-flow  separation,  these  approximations  are  considered  valid. 

b.  Generalized  Axisymmetric  Equations 

The  thin- layer  generalized-axisymmetric  equations  are  obtained 
from  the  three-dimensional  equations  by  making  use  of  two  restrictions: 
(1)  all  body  geometries  are  of  an  axisymmetric  type;  (2)  the  state 
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variables  and  the  contravariant  velocities  do  not  vary  in  the 

circumferential  direction.  In  what  follows,  the  term  of  Eq.  (1) 

shall  be  reduced  to  the  source  term  of  the  generalized  axisymmetric 
equations. 
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A  sketch  of  a  typical  axisymmetric  body  is  shown  in  Figure  la. 

In  order  to  determine  the  circumfer¬ 
ential  variation  of  typical  flow  and 
geometric  parameters,  we  first 
establish  correspondence  between  the 
inertial  Cartesian  coordinates  (x,y,z) 
(to  which  the  dependent  variables  are 
referenced),  the  natural  inertial 
cylindrical  coordinates  (x,4>,R),  and 
the  transformed  variables  (C»n,C). 

The  choice  of  the  independent 
variables  £,  n,  C  is  restricted,  as 
shown  in  Figure  lc,  insofar  as  n  must 
vary  as  <J> ,  i.e.,  <j>  =  Cn  (where  C  is  a 
constant).  From  the  views  shown  in 
Figure  1,  the  relationship  between  the 
coordinate  systems  are  observed  to  be 
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P1GURE  1  •  AXISYWETRIC  30DY  AND  COORDINATE  SYSTEM 


4>  =  Cn 
x  =  xU,c,t) 
y  =  R(£,c  ,r)sin<t> 
z  =  R(£,c,x)cos<f> 


(4) 


where  <)>  =  <J>  (t )  and  the  Cartesian  and  cylindrical  coordinates  are 
related  in  the  usual  way.  Note  that  x  and  R  are  general  functions  of 
only  Z, ,  £ ,  and  x  . 

Evaluating  the  metric  terms  given  the  above  assumptions  and 
substituting  in  Equation. (1)  the  resulting  thin- layer  generalized 
axisymmetric  equations ^5'  can  be  written  as 

6  q  +  6_E  +  3  G  +  H  =  Re"1  3  S  (5) 

x  5  Z  £ 


where 
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H  =  J 


0 

0 


PV[RC(U-Ct)  +  RJW-^)] 

-PVR*n(V-nt)  -  p/(R*n) 

0 


(6) 


is  the  resultant  source  term  which  has  replaced  3^F  of  Equation  (1) 

Equation  (5)  contains  only  two  spatial  derivatives  but  does 
retain  all  three  momentum  equations  thus  allowing  a  degree  of  generalit 
over  the  standard  axisymmetric  equations.  In  particular,  the 
circumferential  velocity  is  not  assumed  to  be  zero  allowing  then 
computations  for  spinning  projectiles  or  swirl  flow  to  be  accomplished. 


The  numberical  algorithm  used  for  both  equations  (1)  and  (5) 
is  a  fully  implicit,  approximately  factored  finite  difference  scheme  as 
analyzed  by  Beam  and  Warming18  .  The  algorithm  for  Equation  (1) 
written  in  the  delta  form  is 


(I  +  h6?An  -  (I  +  h6^Bn  -  e  T J ” 1 V^V^J ) x 

(I  +  h6  Cn  -  hRe1 6  J_1MnJ  e  J_1V  V  J)(qn+1  -  q") 

S  (7) 

=  -At (6  En  *  i  Fn  *  5  C"  -  Re'1 6  Sn) 

s  n  s  s 

-  eEJ_1[(V?Vc)2  +  (V^)2  ♦  (VcVr)2].Jqn 


where  the  6's  are  central  difference  operators,  A  and  V  are  forward  or 
backward  difference  operators;  h  =  At  corresponds  to  Euler  implicit 
first-order  and  h  =  At/2  to  trapezoidal  second-order  time  accuracy. 

The  scheme  can  be  first  or  second  order  accurate  in  time  and  second  or 
fourth  order  accurate  in  space.  A  similar  algorithm  can  be  written  for 
Equation  (5)  as 
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(I  +  hS  An  -  e  J_1V  V  J)(I  +  h6  Cn  -  efJ_1V  V  J 

s  A  s  ^  C  A  C  C 

-  hRe"1 6  J_1MnJ)  X  (qn+1  -  q11)  =  -&t(6  En  +  6  Gn  (8) 

Q  s  c 

-Re_16  Sn)-AtUn  -  e  )2]Jqn 

Notice  that  the  second  factored  terms  of  Equation  (7)  has  been  reduced 

to  the  term  AtHn  and  appears  on  the  right  hand  side  of  Equation  (8). 

The  details  of  the  numerical  method,  algorithm  and  boundary  conditions 
for  each  formulation  can  be  found  in  References  (4)  and  (5)  respectively. 

III.  MODEL  GEOMETRY  AND  EXPERIMENTAL  MEASUREMENTS 

The  ogive- cylinder  boattail  shape  used  for  this  study  closely 
resembles  a  modern  low-drag  artillery  projectile  and  can  be  seen  in  the 
pressure  coefficient  plot  of  Figure  2.  The  model  has  a  three  caliber 
ogive,  a  two  caliber  cylinder,  and  a  one  caliber  7-degree  boattail; 
the  model  length  is  343mm  and  the  diameter  is  57.2mm. 

Several  wind  tunnel  experiments  have  been  conducted  for  this 
model  geometry  in  order  to  obtain  data  for  comparison  to  numerical 
computations.  Experimental  data  used  for  comparison  in  this  study  are 
boundary  layer  profile  measurements  obtained  at  the  Naval  Surface 
Weapons  Center  (NSWC),  White  Oak  Laboratory  and  surface  pressure 
measurements  obtained  at  the  NASA  Langley  Research  Center. 

Boundary  layer  data  were  obtained  in  the  NSWC  Tunnel  No.  2 
which  has  an  open  jet  test  section  with  a  nozzle  exit  size  of 
40. b  x  40.6  cm.  Data  were  acquired  at  Mach  =  0.908  with  a  supply 
pressure  of  one  atmosphere  and  a  supply  temperature  of  approximately 
320°k;  these  conditions  give  a  Reynolds  number  of  4.5  x  10P  based  on 
model  length.  The  NSWC  l.aser  Velocimeter  used  to  measure  the  boundary 
layer  velocities  is  commonly  referred  to  as  a  forward  scatter, 
differential  Doppler  or  "fringe"  system  and  is  described  in  Reference 
(9).  Previous  boundary  layer  measurements  at  Mach  =  3.0  obtained  by 
both  laser  velocimeter  and  by  impact  probe  measurements  showed  good 
agreement.  The  present  Mach  =  0.908  profiles,  however  were  not  as 
consistently  smooth  and  uniform  as  the  Mach  3  profiles  described  in 
Reference  (10). 
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Surface  pressure  measurements  were  obtained  in  the  NASA 
Langley  8  ft  Pressure  Tunnel  with  supply  pressures  and  temperatures  of 
one  atmosphere  and  320° K  respectively.  The  model  was  instrumented 
with  pressure  taps  at  15  longitudinal  stations.  Data  were  acquired 
at  Mach  numbers  of  0.91  to  1.20,  angles  of  attack  from  0  to  10  degrees, 
and  for  circumferencial  positions  by  rolling  the  model  in  22.5  degree 
increments. 


IV.  COMPARISONS  BETWEEN  COMPUTATION  AND  EXPERIMENT 


Computations  have  been  obtained  on  the  same  secant-ogive- 
cylinder  boattail  (SOCBT)  shape  used  for  the  experiment.  The  3-D 
Navier-Stokes  equations  were  used  for  a  >  0°  and  the  generalized 
axisymmetric  formulation  for  a  =  0°.  Computational  results  have  been 
obtained  for  a  ring  airfoil  shape  which  show  an  interesting  shock 
pattern  being  developed  as  a  function  of  Mach  number. 

a.  Secant-Ogive-Cylinder- Boattail ,  a  =  0° 


A  computational  finite  difference  mesh  of  78  longitudinal 
points  by  50  normal  points  were  used  for  these  calculations.  For  the 
test  condition  of  Mach  =  0.908  the  surface  pressure 
C  is  shown  as  a  function  of  axial  position  in  Figure  2. 

P  The  computational  results, 

indicated  by  the  solid  line  are 
shown  to  be  in  good  agreement  with 


experimental 

coefficient. 


the  experimental  results, 
expansion  (decrease  in  C^) 

subsequent  shock  (increase 
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can  be  seen  occurring  near  the 
nose-cylinder  and  cylinder-boattail 
junctions.  An  accurate  computation 
of  the  pressure  distribution  is 
important  since  it  is  the 
integration  of  surface  pressures 
which  are  primarily  used  to 
determine  the  aerodynamic 
coefficients . 


Increasing  the  Mach  number  to  M  =  0.96  (Figure  31  shows  the 
movement  of  shock  position,  which  is  typical  of  transonic  flow. 

Again  the  experiment  and  computations  show  good  agreement  with  the 
exception  of  the  boattail  region.  This  discrepancy  is  attributed  to 
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FIGURE  4.  mat  CONTOURS  M  M  -  .96  TOR  SOCBT,  a  -  0.0 

inadequate  grid  resolution  in  the  vicinity  of  the  shock.  A  more 
detailed  look  at  the  computed  shock  position  can  be  seen  in  Figure  4 
(above)  where  Mach  contours  have  been  plotted.  The  coalesence  of  the 
Mach  lines  represent  the  position  of  the  shock  and  regions  of  subsonic 
and  supersonic  flow  are  identified. 

The  accurate  determination  of  aerodynamic  coefficients  is 
extremely  important  throughout  the  transonic  regime  since  the 
magnitude  of  the  coefficients  can  change  by  as  much  as  100%  in  this 
area.  A  series  of  computations  were  obtained  from  M  =  0.8  through 
M  =  1.1  and  the  surface  pressures  were  integrated  to  determine  the 
aerodynamic  wave  drag.  The  results  are  shown  in  Figure  5  together  with 

the  experimental  data.  Excellent 
agreement  is  shown  for  both  the 
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The  dynamic  stability  of  shell  is  one  area  of  concern  when 
designing  new  shell  or  modifying  existing  ones.  The  Magnus  moment, 
which  affects  the  dynamic  stability,  is  a  viscous  phenomena.  Therefore 
an  accurate  representation  of  the  viscous  portion  of  the  flow  field  is 
crucial  to  computing  the  Magnus  moment.  As  an  initial  attempt  to  look 
at  the  boundary  layer  in  transonic  flow,  with  its  associated  shock 
interaction,  computations  were  performed  at  M  =  0.908,  a  =  0°  and 
comparison  of  computational  velocity  profiles  with  experimental  profiles 
were  made;  an  example  of  this  comparison  is  shown  in  Figure  6.  The 

Navier-Stokes  result  is  the 
solid  line  and  the  circles  are 
the  experimental  results. 

Although  the  computational  and 
experimental  velocity  profiles 
are  not  in  particularly  good 
agreement,  the  shape  of  the 
computed  profiles  is  seen  to  be 
characteristic  of  a  turbulent 
boundary  layer.  A  power  law 
velocity  profile,  characteristic 
of  turbulent  flow,  is  included 
for  comparison. 


b.  Secant-Ogive  Cylinder  Boattail,  a  =  2° 

A  new  finite  difference  mesh  consisting  of  60  longitudinal 
points,  28  normal  points  and  20  points  in  the  circumferential 
direction  was  designed  for  computations  at  angle  of  attack.  Clustering 
of  the  longitudinal  points  was  maintained  in  the  vicinity  of  the 
expansion  similar  to  the  a  =  0°  cases.  The  computed  leeward  (<J>  =  0°) 
surface  pressure  coefficients,  as  a  function  of  longitudinal  position 
are  shown  in  Figure  7  compared  to  the  experimental  results.  The 
computed  results  are  shown  to  follow  the  same  trend  as  the  experimental 
data.  The  agreement  in  the  vicinity  of  both  expansions  is  quite  good 
but  falls  off  on  the  cylinder  and  boattail.  A  major  problem  in  running 
3-D  computations  is  the  large  number  of  grid  points  required  for 
adequate  resolution.  Increasing  the  number  of  longitudinal  points,  similar 
to  the  amount  used  for  a  =  0°  (78  points)  should  result  in  better 
agreement . 
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FIGURE  7.  COMPARISON  OF  EJPCRIWfTA’.  A  NO  COMPUTATIONAL  SURFACE 
PRESSURE  COEFFICIENTS  FOR  SOCHT.  »  •  2« 
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c.  Hollow  Projectile,  a  =  0° 


Of  current  interest  in  shell  design  is  the  utilization  of 
hollow  projectiles  which  have  the  characteristic  of  "flat"  trajectories. 

A  shape  of  this  type,  known  as  the  ring  airfoil,  has  been  type 
classified  and  is  currently  used  as  an  anti  riot  device.  In  order  to 
demonstrate  the  general  geometry  capabilities  of  the  Navier-Stokes  codes, 
computations  were  performed  for  a  ring  airfoil  shape  at  a  =  0° .  A 
cross  section  of  the  actual  shape  is  shown  at  the  bottom  of  the  C  plot 


of  Figure  8a.  Inviscid  results  are  presented  in  Figures  8a,  b,  c  and  d 
for  M  =  0.4,  0.7,  0.8  and  0.9.  In  all  cases  the  pressure  distribution 
is  plotted  for  the  internal  and  external  surfaces  using  a  solid  line 
and  dashed  line  respectively. 
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The  critical  pressure  coefficient,  C^*,  is  the  value  of  the 

pressure  coefficient  at  sonic  velocity.  The  flow  velocity  in  regions 
with  pressure  greater  than  C  *  is  subsonic  and  in  regions  where  the 


pressure  is  less  than  C^*  the  flow  is  supersonic.  For  M  =  0.4 

(Figure  8a),  all  values  of  C  are  greater  than  C^*  indicating  the  flow 

over  both  the  internal  and  external  surfaces  is  subsonic.  However,  as 
the  Mach  number  is  increased  to  0.7,  the  upper  surface  is  shown  to 
develop  supersonic  flow  and  a  shock  wave  while  the  internal  flow 
remains  subsonic.  Increasing  the  Mach  number  still  further  to  M  =  0.8, 
the  shock  waves  are  now  seen  to  exist  on  both  the  external  and  internal 
surfaces.  A  final  solution  at  M  =  0.9  shows  that  the  shock  wave  has 
moved  to  the  trailing  edge  indicating  supersonic  flow  over  most  of  the 
internal  and  external  surfaces. 


V.  SUMMARY 


Implicit  finite  difference  methods  have  been  used  to  solve 
the  thin- layer  Navier- Stokes  equations.  Both  the  three-dimensional  and 
generalized  axisymmetric  equations  have  been  presented  and  solutions  of 
the  flow  field  about  projectile  shapes  have  been  obtained. 

The  computed  surface  pressure  coefficients  on  the  ogivc- 
cylinder-boattail  projectile  were  found  to  be  in  excellent  agreement 
with  the  experimental  data  for  a  =  0° ,  The  generalized  axisymmetric 
equations  used  for  these  computations  are  similar  to  the  2-D  equations 
and  thus  have  no  severe  limitation  on  the  number  of  mesh  points 
required  for  a  good  solution. 

Computations  for  the  same  projectile  shape  at  a  =  2°,  using 
the  3-D  equations,  show  good  agreement  with  the  experimental  data. 
Computational  experimentation  indicate  that  improved  agreement  could  be 
obtained  with  increased  grid  resolution. 

Computations  of  the  viscous  boundary  layer  indicate  the 
correct  trend  for  a  turbulent  velocity  profile  which  is  critical  to 
the  determination  of  Magnus  moment.  Additional  experimental  velocity 
profiles  are  required  to  fully  access  the  accuracy  of  the  Navier-Stokcs 
computations  in  the  viscous  region. 

The  general  geometry  capability  of  the  numerical  scheme  was 
demonstrated  by  solving  the  inviscid  flow  field  about  a  ring  airfoil 
projectile.  The  ability  to  compute  in  regions  of  pure  subsonic  flow 
and  mixed  subsonic/supersonic  flow  has  been  demonstrated. 
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